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In the rapidly emerging field of nanotechnology, as well as in biology where chemical reaction 
phenomena take place in systems with characteristic length scales ranging from micrometer to the 
nanometer range, understanding of chemical kinetics in restricted geometries is of increasing interest. 
In particular, there is a need to develop more accurate theoretical methods. We used many-particle- 
density-function formalism (originally developed to study infinite systems) in its simplest form (pair 
approach) to study two-species A-\- B ^ reaction-diffusion model in a finite volume. For simplicity 
reasons, it is assumed that geometry of the system is one-dimensional (Id) and closed into the ring 
to avoid boundary effects. The two types of initial conditions are studied with (i) equal initial 
number of A and B particles A'^o,^ = No,b and (ii) initial number of particles is only equal in 
average {No, a) ~ {No,b}- In both cases it was assumed that in the initial state the particles are 
well mixed. It is found that particle concentration decays exponentially for both types of initial 
conditions. In the case of the type (ii) initial condition, the results of the pair-like analytical model 
agrees qualitatively with computer experiment (Monte Carlo simulation), while less agreement was 
obtained for the type (i) initial condition, and the reasons for such behavior are discussed. 
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I. INTRODUCTION 



Classically, biochemical reaction kinetics is extrapo- 
lated from measurements in dilute solutions and fitted 
into the cellular reaction environment, and several flaws 
in this approach have been pointed out.^ The main mo- 
tivation of our work is to improve the understanding 
of diffusion-controlled reactions in topologically complex 
nanoscale environments represented in biological cells. In 
this study attempt is made to develop better theoretical 
methods which could describe diffusion-controlled reac- 
tions with boundaries. To achieve this goal, a particu- 
lar way of doing calculation, the many-particle-density- 
function (MPDF) approach, ^"'^ is modified to account for 
presence of boundaries. 

The theoretical findings of this study are relevant for 
experimental work done in refs. 6-8. Even if we focus on 
biochemical reaction kinetics, the results should have an 
equal bearing on nanotechnological applications such as 



nanofiuidics^ or molecular electronics^". Both are likely 
to be strongly dependent on reaction-diffusion behaviors 
of molecules (nanofluidics) or electrons and holes (molec- 
ular electronics) in restricted nanoscale geometries. 

Most of the studies on diffusion-controlled reactions 
have been performed for infinite systems without bound- 
aries and a variety of methods have been developed to 
do the analysis. The methods range from mean field 
treatments towards more exact approaches which em- 
ploy quantum spin-chains^^, field theory^^'^'^, or MPDF 
formalism^"^ . The references 3, 14 are an excellent re- 
view on the subject. The opposite case when reactions 
take place in restricted geometries with reactants con- 
fined into finite size, and eventually squeezed into very 
small volumes, is less understood. There is, however, 
some pioneering work in this area.^^^^^ In here, focus is 
on testing performance of MPDF on diffusion-controlled 
reactions in finite volumes. 

The infinite diffusion-controlled systems posses quite 
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remarkable properties. When dimensionality of the sys- 
tem d is lower than some critical dimension dc (e.g. for 
A+A reaction dc 2^^, and for A + B dc = 4^") a new 
non-trivial sort of kinetics sets in. Taking ^ + 5^0 
as an example: Classical chemical kinetic rate equation 
for this reaction, with initial densities equal and homoge- 
neous, is given by h{t) = —\n{t)'^ where dot over symbol 
denotes time derivative. For large time t this equation 
would predict density decay in the form of power law 
n{t) « At~". The amplitude of decay equals A = 1/X 
and decay exponent is given by a = 1. In reality, a = 1 
holds only for a sufficiently high dimensionality of the 
system when d > 4, while for d < 4 one has a ~ d/A 
and A — const ^/noD^'^/^, and const is just a numerical 
factor. Please note that for d < 4 decay amplitude A 
does not depend on the reaction rate A and exhibits lack 
of dependence on chemical details (universality). 

The kinetics of the type described above is commonly 
referred to as anomalous or fluctuation- dominated. The 
term "anomalous" points to the fact that mean field (or 
classical rate equations) fail to describe such systems. 
The phrase "fluctuation-dominated" emphasizes impor- 
tance of fluctuations in particle densities. Once the reac- 
tion creates a hole in the particle concentration, diffusion 
is very slow in restoring the homogeneous particle den- 
sity. This has to do with recurrence of random walks. 
For d < 2 probability that the random walker will re- 
turn to the same site after arbitrary number of steps is 
equal to one. Random walkers tend to wander around 
their initial position, and particles do not mix that well. 
Rule of thumb is that for lower dimensions the kinetics 
gets more anomalous. Role of dimensionality is well un- 
derstood both for integer and fractal (non-integer) like 
dimensions. On the other hand, much less is known 
what happens when one shrinks the system size, which 
is studied here. 

To impact some progress in understanding reactions 
in restricted geometry we analyze performance of MPDF 
approach and modify it to account for flnite reaction vol- 
ume. To test such method of calculation the A-l-B model 
is used as a study case. The A+B model is natural choice 
for such a task. This model has been intensively studied 
for inflnite system sizes. ^^^"'^^ ^^ It was found that the 
A-l-B reaction has the remarkable property that domains 
rich with A or B particles are formed as time goes on. 
Once domains are formed the reactions happen only at 
domain boundaries which leads to already mentioned de- 
cay exponent a = d/4. It is interesting to study in what 
way the dynamics of the system (kinetics) changes as one 
reduces the volume available to reaction, in particular 
whether domain-like structure survives. For simplicity 
reasons, the one-dimensional (Id) case is studied. In the 
calculation that follows there is nothing special about Id 
and present analysis can be easily extended to the two- 
or three- dimensional cases. 

The A-l-B model in restricted geometry has been stud- 



ied before with the assumption that one type of reactant 
is attached at the center of a small volume, and it was 
further assumed that one type of particle is in large ex- 
cess. The more realistic problem where all particles 
are allowed to move is much harder to solve, and the 
goal of our study is to describe such a situation. Also, 
in here, the focus is on the case when the initial number 
of reactants is the same, or roughly the same. Naturally, 
the shape of the reaction container might be important 
but this issue is not addressed at the moment. To avoid 
boundaries completely, our Id system will be closed into 
the ring. 

The paper is organized as follows. In section II the 
model is developed, i. e. detailed account is given of how 
particles move and react. Lattice model is used due to 
its conceptual simplicity. In section III equations of mo- 
tion are derived using MPDF formalism in its simplest 
form (pair approach). In section IV equations of mo- 
tion are solved analytically and it is shown how multi- 
exponential decay emerges. The results of computer ex- 
periments (Monte Carlo Simulations) are given in section 
V followed by a comparison between theory and Monte 
Carlo simulations in section VI. We conclude by analysis 
of strengths and weaknesses of the pair-approach applied 
to a reactions in restricted geometry in section VII. In 
appendix A details are given how to calculate effective 
reaction rate k{t) which determines density decay. In 
appendix B the algorithm used to do Monte Carlo simu- 
lations is discussed in detail. 



II. THE LATTICE MODEL 

To test any theory one inevitably needs a model which 
serves as study case. The model used here is defined 
as follows. The two species, A and B, move on a Id 
lattice performing random jumps with rates (diffusion 
constants) Da and Db respectively. It is assumed that 
Da — Db- Position of lattice sites is given by Xi = ih 
with i = 0, 1, 2, 3, • • • , M and h denotes lattice spacing. 
Sometimes, x and y will be used instead of Xi. Periodic 
boundary conditions are assumed and sites « = and 
i = M are defined to be equivalent. There are M lattice 
sites in total and L = Mh. By using periodic bound- 
ary conditions it is possible to work with a system of 
finite size and yet keep the spatial translational invari- 
ance. This greatly facilitates the analytical treatment of 
the problem. 

It is assumed that the reaction probability (per unit 
time) for particle A a.t x and S at y is given hy a{x — y). 
For a{x — y) simplest possible form is used 

(j{x-y)=(7(i6{a-\x-y\). (1) 

where 6(x) = for x < and 9{x) = 1 for a; > 0. In 
this way two important aspects of chemical reactions are 
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embedded, a corresponds to the effective range of reac- 
tion and (Jo is its strength. One could also say that each 
particle carries a ring of radius a/2 and when two rings 
overlap the particles can react. In this sense a/2 could be 
thought of as the size (radius) of particles. For simplic- 
ity reasons it is assumed that the reaction products do 
neither influence reactants, nor the A-|-B reaction. Also, 
exclusion or steric effects are not taken into account, i.e. 
particles are allowed to "enter" into each other (please see 
Fig. 1) and react with same probability independently 
from which direction they approach each other. 

The model has the useful property that if a is thought 
of as the size of reactants, then by varying a several in- 
teresting situations can be studied. For example, when 
o is on the order of the system size L, one can think 
of situations of extreme crowding. On the other hand 
when a <^ L reactants appear as point-like objects. In 
Fig. 1 we offer a schematic way how to think about these 
situations. The model presented above is solved analyti- 
cally and numerically by a Monte Carlo simulation in the 
following sections. 



III. EQUATIONS OF MOTION IN PAIR 
APPROXIMATION 

To solve the A-l-B reaction-diffusion model in a 
restricted geometry we use a many-particle-density- 
function formalism (MPDF), since it was already used 
to describe asymptotics of the same reaction in an infi- 
nite volume. We modify the formalism and apply it 
to the case of a restricted geometry. In the following the 
formulation presented in ref 2 will be closely followed. 
On the way, the changes made to the original formalism 
will be discussed. 

The dynamics of the system, as defined in previous 
section, is stochastic and governed by Master Equation 
which describes time evolution of conflgurational proba- 
bilities of the system P{c, t), 

He, t) = Yl [Wc'^cP{c', t) - Wc^c'P{c, t)] (2) 

c' 

where c is short notation for occupancy of lattice sites 
and Wc^rj arc transition probabilities which can easily 
be deduced from the previous description of the model. 
Here and throughout the paper dot over symbol denotes 
time derivative. 

The quantities of interest are particle densities pa{x, t) 
and pBix,t) and they can be calculated from P{c,t) (at 
least in principle). Since system is closed into a ring 
translational invariance holds and concentrations cease 
to be position dependent which leads to pa{x, t) = nA(t) 
and pB{x,t) = nsit). Following recipe in ref. 2 gives 
following equations for ua and ris. 



fL/2 

fiA{t) = -nA{t)nB{t) dxa{x)Y{x,t) (3) 

J-L/2 
fL/2 

hB{t) = —nA{t)nB{t) I dxa{x)Y{x,t) (4) 

J-L/2 

The Y{x,t) denotes correlation function for AB pairs. 
Absence of correlations is signaled hyY{x,t) = 1. Please 
note that in this work system size is finite which enters 
through finite integration domain in the integrals above 
(it might appear as minor technical detail but this fact is 
very important). Also, it is assumed that reversal sym- 
metry holds, i.e. Y{x,t) = Y{—x,t). 

Again, following ref. 2 one can derive equation for 
Y{x, t) which is given by 

Y{x, t) = {Da + Db)Y"{x, t) - a{x)Y{x, t) 
rL/2 

-nBY{x, t) / dya{y)Y{y, t) [Xb{x -y,t)~ 1] 

J-L/2 
rL/2 

-nAY{x, t) / dya{y)Y{y, t) [Xa{x - y,t) - 1] (5) 

J-L/2 

where prime denotes spatial derivative, and Xa{x, t) and 
XB{x,t) correlation functions for AA and BB pairs re- 
spectively. The Xa{x, t) and Xb{x, t) obey similar equa- 
tions which are not given here to save the space. 

The equations (3)- (5) are derived under assumption 
of Kirkwood superposition approximation, which is a 
technical way of saying that dynamics is governed by 
pair effects. Naturally, assumption of the dominance of 
pairs effects is an approximation. It might or might not 
work, and the goal of present study is to test this. In 
the following, to make analytic treatment possible, equa- 
tions will be simplified further by setting XA{x,t) and 
XB{x,t) equal to one. This amounts to ignoring corre- 
lations among AA and BB pairs. In ref. 2 it was shown 
that (for infinite reaction volume) such approximation 
is too severe and does not lead to correct decay expo- 
nent a = d/4 (it gives a = d/2). Nevertheless, in here 
we consider such simplification. The validity of such an 
approximation, together with the fact that we are using 
Kirkwood approximation, is tested via computer experi- 
ment later on. 

The form of boundary conditions for Y{x,t) differs 
from the one used in ref. 2. In the case of infinite system 
one takes 

Y{x,t) ^1 , x^oo (6) 

while for finite system with periodic boundary conditions 
another form has to be used 

Y{x + L,t)=Y{x,t) (7) 

It will be shown later that the change from (6) to (7) 
leads to a qualitative change from power law to (multi) 
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exponential behavior for correlation dynamics. The rest 
of the boundary conditions are standard, and are taken 
as in the case of an infinite system size, 



nA{0) = no 
nsiO) = no 
Y{x,0) = 1 



(8) 
(9) 
(10) 



Also, taking L ^ oo should reproduce findings of ref. 2, 
within set of approximations employed here. 



IV. EMERGENCE OF MULTI-EXPONENTIAL 
DENSITY DECAY 



With assumptions Xa = Xb = 1 Eq. (5) reduces to 

Y{x,t) = {Da + DB)Y"{x,t) - a{x)Y{x,t) (11) 

Eq. (11) is solved by using a Laplace transform as shown 
in the appendix A. To simplify the algebra, it is as- 
sumed that o"o is arbitrary large. The exponential behav- 
ior emerges due to the fact that the spectrum of Eq. (11) 
is discrete due to particular nature of the boundary con- 
ditions. The final expression for k{t) reads 



k{t) = kreg{t)+2ad{t) 



(12) 



and details of calculation are given in the appendix A. 
The regular part of k{t) is given by 



L - 2a 



(13) 



1=1 



and Km are constants of multi-exponential decay (eigen- 
values). 



f2m- 1 
\L-2a 



(14) 



The i5-function term in (12) arises from the second term 
on the right hand side of Eq. (A5) when ctq — > oo (please 
see the appendix A). 

Once k{t) is available one can calculate n{t) as 

n{t) = .""^ „ (15) 

^ ' 1 + I{t)no + 2ano ^ ' 

where I{t) = Jq dt'kreg{t') and 

m—l 

The 2ano term in the denominator of (15) comes from the 
d{t) term in Eq. (12). It describes the immediate annihi- 
lation of particles which are within reaction range. When 



(To ^ oo, this happens instantaneously. Thus there is a 
sudden jump in particle concentration. For finite (Tq this 
jump becomes a smooth transition (exponential decay 
with decay exponent proportional to large number cto)- 

The question is whether one can obtain results for an 
infinite system from the Eq. (15) above. This can be done 
using Poisson resummation formula when Dt/{^ — a)'^ <C 
1. The Poisson resummation procedure gives I{t) ^ t^/^ 
which results in the wrong exponent for the density de- 
cay; n{t) ^ instead of correct n{t) ^ t"'^/^. Thus 
we just reconfirm the well known fact that for infinite sys- 
tems, the pair approach predicts too fast decay of par- 
ticles. However, for finite systems, the situation is not 
that clear, it appears to depend on the type of initial 
conditions the real system is subject to. 

In here we consider two types of initial conditions, (i) 
When initially there is an equal number of A and B parti- 
cles, one has n — > as t — > c»; and at the end all particles 
have to annihilate, (ii) One can look at an ensemble of 
similar systems with equal number of A and B particles 
at f = on average, {Nq^a) = (A^o,b)- In such a case, one 
has different asymptotics, {N{t)) — > N{oo) ast^oo. 

Theoretical prediction is that, as time goes to infinity, 
the particle density exponentially approaches the value 
nth{oo); 



nth{oo) = 



no 



1 + noL 



(17) 



The value for nthioo) above can be obtained by sending 
t — > 00 in Eqs. (15) and (16). From (17) one sees that 
asymptotically number of particles is given by 



Nth{^)=No/{l + No) 



(18) 



where Nth{oo) = Lnth{oo). Please note that Nth{oo) 
never approaches zero and settles at a number between 
zero and one. In the case of type (i) initial conditions, 
all particles annihilate and N{oo) = 0. This is clearly 
in contradiction with Eq. (17). However, situation is not 
that hopeless, as will be discussed later. For the type (ii) 
initial condition, for each member in ensemble there is a 
chance that more than one particle will be left, since one 
start dynamics with (random) excess of A or B particles 
at t = 0. Thus, in average, (iV(oo)) will be larger than 0. 
Clearly, pair approach has more chance to describe this 
situation correctly. 

In summary, we find an exponential decay in the long 
time limit which is a pure artifact of the finites of the 
system. There is a clear indication that the quality of 
prediction depends on the type of initial conditions used 
in experiment. Also, the approximations made in deriv- 
ing Eq. (15) are rather severe and in order to check the 
applicability of such a pair-approach Monte Carlo simu- 
lation is used. 
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V. RESULTS OF MONTE CARLO SIMULATIONS 
OF A+B REACTION IN RESTRICTED 
GEOMETRY 



Figures 2,3, and 4 summarize the results of the Monte 
Carlo simulations in. d = 1. The Monte Carlo algorithm 
is described in detail in appendix B. Figure 2 shows a 
simulation for a system with a large initial number of 
particles with a varying reaction range from a nearest- 
neighbor interaction with a/L = 0.0001 towards a longer 
range with a/L = 0.02. Figure 3 shows the case when 
there are initially very few (exactly 10=5A+5B) particles 
present in the reaction volume, also with varying reaction 
ranges from a/L = 0.001 to a/L = 0.2. Thus figures 2 
and 3 give simulation results for type (i) initial condition. 
Figure 4 deals with type (ii) initial condition, when the 
initial number of particles in an ensemble varies with the 
constrain that the sum of A and B particles equals 10. 
(For example, one run could be done with 7A and 3B, 
the other nm with 5 A and 5B, and a third run with 4 A 
and 6B particles, etc.) Figure 1 is a sketch of how to 
think of various situations when a changes from small to 
large values. 

From figure 2 it can be seen that in the case of the near- 
est neighbor reaction range (n = 1) four distinct regimes 
appear and the log-log plot is used to reveal them; (a) 
mean field decay, (b) the plateau region (c) power law de- 
cay and (d) exponential decay at the end. These regimes 
disappear as the reaction range is increased, and eventu- 
ally, for very large a, one only has the exponential regime. 

The mean field regime corresponds to annihilation 
of particles with all reactants being well mixed. This 
leads to depletion of lattice to concentration of the or- 
der n ^ 1/a, thus one particle per reaction range. Then 
diffusion starts to operate and mixes particles. What is 
interesting is that for very large values of o, the plateau 
region starts earlier and lasts longer. Apparently, it takes 
some time before the particles find each other by diffusion 
and start reacting again. 

The power law decay starts after the plateau region. 
There is universality in the power-law regime since all 
curves with different values of a merge into one. This is 

somewhat surprising since a larger a should mean faster 
annihilation, which indeed happens in the mean field 
regime, but yet in the power law stage all curves share 
the same power-law behavior. We speculate that this has 
to do with self organization and build up of correlation. 

The exponential regime is entered after the power law 
regime, when the number of particles in the system be- 
comes small. With the present computer hardware it 
was not possible to resolve this exponential regime bet- 
ter. This is indeed done in figure 3 with a smaller lattice 
size and lower particle number. 

To illuminate this exponential decay at later stages of 
annihilation, we performed simulations with a smaller 



number of particles (10=5AH-5B) on a smaller lattice 
with 10^ sites. Thus we used type (i) initial condition. 
To obtain each curve we followed 1000-3000 realizations 
of dynamics and averaged over such an ensemble. The 
result is shown in figure 3. The upper figure is in log-t 
scale to resolve the small and large t region, respectively. 
The lower figure is in normal-t scale and we use it to 
detect exponential decay (where a straight line indicates 
exponential decay). 

The crossover from mean-field to plateau-like dynam- 
ics can be seen in the upper graph where all curves 
drop down to a plateau value which is a-dcpendent. 
The theoretical prediction for this plateau is n{0~^) = 
n-o/i^ + 2ano). The initial drop in concentration is large 
for large a-values. In the upper figure, it is hard to say 
when the plateau behavior turns into exponential decay. 

The lower graph shows that decay is indeed expo- 
nential since density curves at late times are straight 
lines in the logio(n)-t plot. Thus at the late times 
n ~ exp(— Kii). Also, the decay constant ki is a- 
dependent since slopes are different for various values of 
o, and Ki becomes larger with larger a. Also, it appears 
that there is an upper limit for a at which decay becomes 
infinitely fast. Naturally, this happens when a = L/2 
since none of the particles can escape from each other. 
The qualitative dependence of ki on a just discussed is 
in agreement with theoretical prediction in Eq. (14) with 
m = 1. 

Figure 4 is obtained in a similar way as figure 2. The 
only difference is that figure 4 deals with the type (ii) ini- 
tial condition. For the particular run, when Nq^a Nq^b, 
the final number of particles in the system is not zero. 
For example, when starting from 7A and 3B particles, 
the system will end up in the state of 4 A particles. This 
comes from that fact that the A-l-B reaction conserves 
the particle difference A^a(0 — A^s(i) = const. Curves 
for different values of a saturate at one single value which 
is, independent of a. Clearly, the value of the plateau is 
solely controlled by the excess of particles at t = 0, and 
can be calculated from theory if needed, but result of 
Monte Carlo simulation is equally informative. 



VI. COMPARISON OF COMPUTER 
EXPERIMENT AND THEORY 



Figure 5 shows a comparison of the anal}d;ical treat- 
ment with computer experiment (simulation parameters 
as in figure 3). It can be seen that the pair (Smolu- 
chowskii) approach does not predict that the number of 
particles in the system should approach zero. The rea- 
sons for this are discussed later but have to do with the 
fact that we are looking at highly symmetric situation 
with equal number of A and B particles all the time. To 
enforce such zero asymptotics by hand we use interpola- 
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tion formula 

riintit) = n{t) - n{oo){l - e-^^') (19) 

where Ki is the first dominant large time exponent in ex- 
pression for k{t) (see Eqs. 15 and 16). It can be easily 
seen that the equation above holds exactly for f = and 
t ~ oo. If Eq. (19) is used instead of (15) the agreement 
with simulation improves in the sense that the decay is 
exponential and the qualitatively theoretical exponent is 
roughly the same as the one obtained from simulations. 
More work is in progress to develop improved interpola- 
tion formulas. 

Figure 6 deals with the same type of comparison, but 
with a simulation setup as in figure 4 when the initial 
number of particles is not fixed, just the total num- 
ber (type (ii) initial condition). One can see that the 
agreement between theory and simulation is much bet- 
ter. Clearly, Smoluchowskii theory deals better with the 
type (ii) initial condition represented by the figure 4 than 
by the type (i) represented by the figure 3. Also, in figure 
6, one can see that theory predicts too fast particle anni- 
hilation. This is no surprise since this is what one would 
expect form such pair approach which does not take into 
account formation of domains. (In that respect there is 
similarity with infinite systems, but only for the type (ii) 
initial condition). 



VII. DISCUSSION 

The goal of present work was to impact some under- 
standing of diffusion-controlled reactions in restricted ge- 
ometries, with aim to describe some aspects of chemical 
reactions in biological cells. Two issues have been dealt 
with: 

(1) The particular way of doing calculation was tested, 
the MPDF formalism developed by Kuzovkov and Ko- 
tomin (see ref. 2 for details). To be able to solve equa- 
tions analytically the hierarchy of many-particle-densities 
was truncated at the level of three-particle- density us- 
ing shortened Kirkwood superposition approximation. 
This approximation amounts to assuming that pair ef- 
fects dominate correlations, and present calculation can 
be viewed as a variant of pair approach. 

(2) A two species reaction-diffusion model A + B ^ 
in a restricted geometry was taken as a study case. Two 
types of initial conditions were considered, type (i) where 
the initial number of A and B particles is strictly equal, 
and type (ii) initial conditions where the initial number 
of particles is equal only approximatively. 

Thus the paper is best viewed as a method paper since 
the main goal is to test the strengths and weaknesses of 
the pair approach. To test quality of approximations in- 
volved all results have been compared with the results of 
computer experiment (Monte Carlo simulation). 



From a theoretical point of view, it seems that the pair 
method, being widely used in calculation of bulk proper- 
ties, works with mixed success for the restricted reaction 
diffusion systems, at least the one studied in here. The 
agreement between theoretical calculation and computer 
experiment is qualitative in the case of type (ii) initial 
conditions. In the case of type (i) initial conditions there 
is less agreement, however, situation is not that hopeless. 

In the case of type (i) initial conditions pair approach 

makes error of the order of one particle (please see 
Eq. 18), since it predicts that, in the average, after very 
long time, there will be between zero and one particle in 
the system (though all particles should vanish). When 
initial number of particles is relatively large, the pair ap- 
proach can describe evolution of system for rather long 
time, before the regime is reached where only one parti- 
cle is left in the reaction volume. However, in the case 
of small initial number of particles, there is no such time 
interval, and the mismatch between pair approach and 
simulation has to be addressed more seriously. 

The weakness of pair approach in dealing with type (i) 
initial condition rests on the fact that the truncated set 
of equations for many-particle-densities do not recognize 
any effects which go beyond pair correlations. For ex- 
ample, in the present case, all information related to the 
fact that initially there were 5A and 5B particles in the 
system, and that all particles have to vanish eventually, 
is missing. Work is in progress to pass such type of in- 
formation from higher order particle-density- functions to 
lower order ones. For example, we already have better 
interpolation scheme than the one given in Eq. (19) for 
the case when there are initially three particles in the 
system, but we are trying to understand how to extend 
such analysis to higher numbers. 

Interestingly enough, it seems that, contrary to the sys- 
tems with infinite sizes, setting XA{x,t) = Xsix^t) ~ 1 
is reasonable approximation for a finite system, but we 
have to perform more tests. This could have to do with 
the fact that if the system is too small, there will be no 
time to develop clusters of A and B particles, and setting 
XA{x,t) — 1 and XB{x,t) — 1 might turn to be a good 
approximation after all. Thus, one does not have to turn 
to more complicated methods of calculation if qualitative 
results are needed. Nevertheless, it is highly desirable to 
see what happens as one includes correlations among AA 
and BB pairs. 

To be able to solve equations analytically, we had 
to simplify MPDF considerably down to the level of a 
pair like approach and various calculation schemes con- 
tain pair approach as possible approximation. Perhaps 
the most common form of pair approach is the one sug- 
gested by Smoluchowskii (see e.g. ref. 3 for interesting 
review). The Smoluchowskii approach boils down to so- 
lution of Poisson equation with different boundary con- 
ditions. (Many readers will be familiar with this in the 
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context of heat transfer or quantum mechanical prob- 
lems.) 

However, one has to keep in mind that pair approach is 
an approximation, and it has to be tested to see whether 
it works. (For example, the pair approach can not de- 
scribe A-l-B reaction-diffusion model when system size is 
very large. This has been discussed in ref. 4.) The ad- 
vantage of pair approach, in the form used here, is that 
it is possible to go beyond it in a systematic way. 

Also, from the particular way we have approached pair 
problem, one can see that the difficulties associated with 
Eq. (18) are likely to be much deeper than just the fact 
that we are using pair approach. (This problem clearly 
vanishes when system size is infinite, as nth{oo) goes to 
zero.) Any scheme which focuses on low rank particle- 
density-function will suffer in a similar way, in the case 
of the highly symmetric, e.g. type (i), initial condition. 
One really has to find a way how to incorporate infor- 
mation of higher order correlation functions into lower 
order ones, without calculating higher order correlations 
functions explicitly. This is a pressing issue. 

A few words about the model used. The goal of present 
work is to develop calculation method rather than to de- 
scribe specific chemical system. We used the model which 
wc could solve, within reasonable level of approximation. 
Nevertheless, the question whether present model has 
any relevance for real biological and chemical systems 
needs to be addressed. 

The reaction diffusion model studied here appears to 
be too simple, for two reasons. First, it does not account 
for chemical details which enter only through two pa- 
rameters, diffusion constant D and reaction rate A. For 
example, the exclusion effects and steric effects are not 
contained in it. Also, the influence of product molecule 
is completely ignored. Second, Id character of the model 
might be too restrictive. 

Despite its simplicity, the model used here contains 
basic characteristics of diffusion-controlled reactions (re- 
action times tji are much smaller than the correspond- 
ing diffusion times tf)); particles are moving on the on 
the lattice and react when within reaction range with no 
memory of initially velocity. Previous research reviewed 
in refs. 2, 3 has shown that diffusion-controlled models, 
similar to the one used here, can be used to describe real 
chemical reactions. In particular, the A-l-B model have 
been used to study two reactions in capillary tube quite 
successfully, bromine -|- cyclohexene adduct, and Cu^"*" 
-|- disodium ethyl bis (5- tertrazolylazo) acetate trihydrate 
— > 1:1 complex in water. 

Why can one be sloppy and ignore chemical details to 
some extent? The reason for this is universality. Most of- 
ten, predictions of the reaction-diffusion models (on lat- 
tice) are insensitive to the details of the chemistry in- 
volved. (For example, the decay amplitude for A-l-B re- 
action does not depend on the reaction rate A.) This 



statement is valid provided one deals with very large sys- 
tem sizes. Naturally, this view is not the only one. There 
are other ways to approach diffusion-controlled reactions. 
For more chemical or biological approach to diffusion- 
controlled reactions see refs. 14, 31 and 32 respectively. 

The simplicity of the model is not necessarily such a 
big handicap, until one reaches extremely small sizes. For 
small system sizes, density decay will start depending on 
details. However, there is a large window in system sizes 
between extremely large and extremely small where such 
kind of universality could survive. In here we push the 
model over its borders by studying situation of extreme 
crowding and not accounting for exclusion effects. 

Furthermore, we would like to notify the reader that we 
do refer to the model here as a "toy" model, it relatively 
simple to formulate it. However, this does not mean that 
the model is easy to solve, quite the contrary. In the case 
of infinite system size it has taken a lot of research effort 
to clarify that the decay exponent indeed is d/A. This 
issue was finally settled in ref. 27 which provides a strict 
mathematical proof. 

The model has a potential to describe experiments in 
refs. 6-8 where, for example, the average diameter of the 
reaction container (liposome) is L ~ 1 — 25/xm. The re- 
actants A (enzyme) and B (substrate) are of the size of 
as, as ~ Inm and the typical number of reactants in- 
serted is on the order of AT 1000. Thus, a <^ L holds 
to a very good approximation. In these experiments re- 
actants appear as point-like objects and there is no need 
to give structure to reactants. Also, problems associated 
with Eq. (18) will likely not to case any damage due to 
large number of particles at i = 0. 

Prom the data above, the concentration of particles 

UA.B ^ N/L^ , is easily estimated to be n ^ l/im~'^, and 
the typical distance between particles is (Iab ~ 1/c^/^ w 
L/N^/^ ~ O.l/zrn. Thus, in average, as,E (Iab, and 
particles are very well separated. Therefore, it is reason- 
able to assume that pair effects are the dominant ones. 
This in turn simplifies the theoretical description consid- 
erably. Clearly, there are other scales in the problem and 
the criteria on applicability of the pair approach are more 
subtle in reality.^'^ 

To summarize, it would be extremely interesting to 
have a general method of calculation which could describe 
diffusion-controlled reactions in finite volumes, perhaps 
something on the level of the pair approach. Pair ap- 
proach is attractive since inclusion of chemical details 
such as exclusion or steric effects is possible (see e.g. work 
in ref. 33) which certainly opens a new route towards 
more quantitative results. However, pair approach is an 
approximation, and before burdening pair approach with 
increasing amount of chemical details, one has to test its 
ultimate reach. Present work is an attempt in this direc- 
tion. 
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APPENDIX A: DERIVATION OF THE 
REACTION RATE k{t) 



APPENDIX B: COMPUTER EXPERIMENT VIA 
MONTE CARLO SIMULATIONS 



With initial condition Y{x,0) = 1, the Laplace trans- 
form of (11) becomes 

sY{x, s)-1 = {Da + Db)Y"{x, s) - a{x)Y{x, s) (Al) 

The correlation functions are symmetric around origin, 
i.e. Y{x,t) = Y{—x,t). They are also periodic in L. 
This implies that it is sufficient to focus on positive x 
axis and impose boundary conditions ■^Y{x,s) = for 
a; = and x = L/2. Equation (Al) is an ordinary sec- 
ond order differential equation, which is easily solved by 
solving it in regions < x < a and a < x < L/2 sep- 
arately and then matching solutions at the end. After 
some algebra one obtains 

Y,{x,s)^^(--- 
D \v fi 



ch{xy/il) 



ch{a^) + sh(av^)ch[V^(| - a)] 



-F^ (A2) 



and 



Y2{X,S) = 



D 



ch[{^-x)V^ 



+ 



Di 



+ 



(A3) 



where Y{x,s) = Yi{x,s) for < a; < a and Y{x,s) = 
Y2{x,s) for a < X < L/2 with ij, = [s + (Jq)/D and 
V = s/D. The reaction rate k{s) is given by k{s) = 
2(To /q dxYi{x, s) and equals 



sh(ay^) 
ch(a;y^) 



+ (A4) 



ch(a^) + yfsh(a^)ch[^(f - a)] £>M 

We could not find the inverse Laplace transform of the 

expression above in closed analytic form. However, this 
is possible when ctq — > oo. In such a case one has 



1 



s cth[(L/2-a)^ 



+ 



2crna 



CTO 



O{l/ao) 
(A5) 



The inverse Laplace transform of the approximate ex- 
pression for k{s) can be found by a residuum method. 
The s = is not a branching point nor pole. The only 
poles come from cosh term in denominator which has 
poles at Sm = -TT^{m - \/2fL)/{L/2 - af . This fully 
fixes form of fc(t) in Eq. (12). 



We have chosen the minimal process algorithm for the 
simulations for two reasons. The first reason is that the 
algorithm reproduces the master equation (2).^'* Second, 
our goal is to study a whole range of particle sizes and 
relatively large numbers of particles at the same time. 
Clearly, there are another possibilities to carry out Monte 
Carlo simulation, but the main advantage of the minimal 
process algorithm is that it can be applied for systems 
containing relatively large number of particles. An origi- 
nal algorithm was devised for the situation where h, 
i.e. particles react at the same lattice site or when near- 
est neighbors. We had to modify the original version of 
the algorithm to account for finite reaction range when 
a ^ h. A detailed description of the algorithm is given 
bellow. 

Algorithm: 

(1) Site i is chosen at random. 

(2) If the site is empty go to step (5). 

(3) For a chosen site i, one has to calculate the rate Wi 

for a certain process to occur (diffusion or reaction) . 
Also, one needs a null rate A'^; where nothing hap- 
pens (the so called "null process"). The null rate is 
defined from Wi+Ni = Q, where Q is arbitrary but 
known at each simulation step. Q is chosen in such 
a way that none of the Ni is negative. In practise, 
the case when Q is taken as the largest of Wi works 
best since this leads to the smallest possible values 
for Ni, i.e. chance that nothing is done in course 
of simulation is reduced. (Please note that this re- 
quires that Q is updated as simulation proceeds, 
but can be done in a straight forward manner as 
explained in ref. 34) 

Wi = Di + Ri accounts for possibilities that a par- 
ticle at the site diffuses to the neighboring site with 
rate D^, or reacts with a particle in some other site 
with rate Ri = X^jefj; '^i^ij)- denotes set of 
sites which are within reaction range of the site i. 
The calculation of Ri is by far the most costly step 
when a is large. In that case, a large region has to 
be searched in order to find all particles within f7j . 
This step costs Msearch ~ {a/hY computational 
steps if the sites are checked one by one. The cost 
can be reduced further by introducing a list which 
specifies which sites that contain particles within 
the reaction range of the particle at site i. In that 
case one has to update the list for each diffusion 
step made. The best algorithm we have so far up- 
dates the list in roughly Mgerach ^ {a/hY~^ steps. 
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(4) Once the rates for the site i have been calculated 
one can use them to evaluate probabilities for spe- 
cific process p^^^ — Di/Q, p^^^ — Ri/Q and 
pinuii) _ jv/Q. Once the probabilities are calcu- 
lated a certain process is chosen by linear selection 
algorithm. First one decides if diffusion, reaction 
or nothing is going to happen. If diffusion is to 
happen then the particle is moved to one of the 
randomly chosen 2 x d nearest neighbors. If re- 
action was chosen, then one of the sites containing 
particles in reaction range is chosen at random, e.g. 
at site j, and pair of particles from site i and j are 
annihilated. 

(5) Time is updated according to the formula t 

t + At where At = 1/LQ where L was specified 
before and Q is the maximum rate at the present 
step. 

(6) Move back to (1) unless some criteria to stop is 
invoked. 

Applying the same type of reasoning as in ref. 34 one 
can see that the algorithm proposed here reproduces the 
behavior described by the master equation (2). As time 
of the simulation progresses we monitor the number of 
particles and calculate all the statistics. 

As the original minimal process algorithm, the present 
simulation method is not that efficient at the later stages 
of dynamics when the lattice becomes sparse. The quan- 
tity that governs computational cost of this method is the 
number of Monte Carlo steps needed to see some change 
in the number of particles. We describe it by the number 
of Monte Carlo steps needed to annihilate the last pair 
of particles. 

To make such estimate, it is best to move to the 
reference frame of one of those particles. Then one 
particle is fixed and another one is trying to find it. 
The number of diffusion steps that the moving particle 
needs to find the one who sits still is roughly given by 
M(diff) ~ L'^/a (here and in the following it is implic- 
itly assumed that every length variable is measured in 
units of lattice spacing h). Each diffusion step bears 
M(step/diff) computational steps which gives the to- 
tal number of steps to annihilate the pair of particles 
equal to M(tot) = M(diff)M(step/diff). The number of 
Monte Carlo steps per one diffusion step is roughly 1, 
M(step/diff) ~ 1. However, calculation of Ri requires 
updating the internal list which costs Mg^arch ~ a'*"^ 
search steps whenever the particle is moved. Thus, the 
true number of computational steps per diffusion step 
is given by M(step/diff) ~ M search ~ a''"^- Finally, 
one gets an estimate for the number of computational 
steps needed to annihilate the last pair of particles as 
M(tot) ~ L'^a'^-'^. 

The algorithm has an interesting property that for 
d = 1 there is a reduction in the computational cost 



when comparing large and small a cases. For larger a 
the algorithm works more efhciently. For d = 2 the com- 
putational cost does not depend on a. Simulating a large 
a situation for d = 3 is more costly. One could avoid 
this growing cost problem at d = 3 by browsing through 
particles instead of searching for sites when calculating 
Ri. This is clearly the preferred option when the number 
of particles in the system is not that large. 



1 Luby-Phelps, K. Int. Rev. Cytol. 2000, 192, 189. 
^ Kotomin E.; Kuzovkov, V. Rep. Prog. Phys. 1992, 55, 
2079. 

^ Kotomin E.; Kuzovkov V. in Comprehensive Chemical Ki- 
netics, vol. 34, "Modern aspects of diffusion-controlled re- 
actions", R.G.Compton and G. Hancock Editors, (Elsevier, 
1996) 

* Kuzovkov V.N. ; Kotomin, E.A. Chem. Phys. 1983, 81, 335. 
^ Kotomin, E.; Kuzovkov, V.; Frank, W.; Seeger, A. J. Phys. 

A 1994, 27, 1453. 
® Chiu, D.T.; Wilson, C; Ryttsen, P.; Stromberg, A.; Karls- 

son. A.; Nordholm, S.; Hsiao, A.; Gaggar, A.; Garzia- 

Lopez, R.; Moscho, A.; Orwar, O.; Zare, R.N. Science 

1999, 283, 1892. 
'' Chiu, D.T.; Wilson, C.; Karlsson, A.; Danielsson, A.; 

Lundqvist, A.; Stromberg, A.; Ryttsen, F.; Davidson, M.; 

Nordholm, S.; Orwar, O.; Zare, R.N. Chemical Physics 

1999, 247, 133. 
® Karlsson, A.; Karlsson, M.; Karlsson, R.; Cans, A-S.; 

Stromberg, A.; Ryttsen, F.; Orwar, O. Nature 2001, 409, 

150. 

® Karlsson, R.; Karlsson, A.; Karlsson, M.; Cans, A-S.; 
Voinova, M.; Bergenholtz, J.; Ewing, A.G.; Akerman, B.; 
Orwar, O. Langmuir 2002, 18, 4186. 

Park, H.; Park, J.; Lim, A.K.L; Anderson, E.H.; Alivisatos, 
A.P.; McEuen, P.L. Nature 2000, 407, 57. 
" Alcaraz, F.C.; Droz, M.; Henkel, M.; Rittenberg, V. Annals 
of Physics 1994, 230, 250-302. 

Konkoli, Z.; Johannesson, H.; Lee B.P. Phys. Rev. E 1999, 
59, R3787. 

" Konkoli, Z.; .Johannesson, H. Phys. Rev. £; 2000, 62, 3276. 
Comprehensive Chemical Kinetics, Vol. 25, "Diffusion- 
limited reactions", C.H. Bamford, C.F.H. Tipper and R.G. 
Compton Editors, (Elsevier, 1985). 

Khairutdinov, R.F.; Serpone, N. Prog. React. Kinetics 
1996, 21, 1-68. 

Khairutdinov, R.F.; Burshtein, K.Ya.; Serpone, N. J. Pho- 
tochemistry and Photobiology A 1996, 98, 1. 

" Tachiya, M. Chem. Phys. Lett. 1980, 69, 605. 

Photochemistry in Organized and Constrained Media, 
edited by Ramaunirthy, V.; VCH Publishers (1991). 

" Lee, B.P. J. Phys. A 1994, 27, 2633. 

20 Lee, B.P.; Cardy, J. Stat. Phys. 1995, 80, 971. 

2^ Havlin, S.; Ben-Avraham, D. Adv. in Physics 1987, 36, 
695. 



9 



Burlatskii, S.F.; Ovchinnikov, A. A. Russ. J. Phys. Chem. 
1978, 52, 1635. 
2^ Ovchinnikov, A. A.; Zeldovich, Ya.B. Chem. Phys. 1978, 
28, 215. 

2* Toussaint, D.; Wilczek, F. J. Chem. Phys. 1983, 78, 2642. 
"^^ Burlatskii, S.F.; Ovchinnikov, A.A.; Pronin, K.A. JETP 
1987, 92, 625. 

26 Gutin, A.M.; Mikhailov, A.S.; Yashin, V.V. JETP 1987, 
92, 941. 

"^"^ Bramson, M.; Lebowitz, J.L. Phys. Rev. Lett. 1988, 61, 
2397. 

Oerding, K. J. Phys. A 1996, 29, 7051. 
Mattis, D.C.; Glasser, M.L. Rev. Mod. Phys. 1998, 70, 979. 
^° Koo, Y.E.L.; Kopelman, R. J. 5«a*. P%s. 1991, 65, 893- 
918 

Calef, D.F.; Deutch, M. Ann. Rev. Phys. Chem. 1983, 34, 
493. 

^2 Berg, H.C.; Purcell, E.M. Biophys. J. 1977, 20, 193. 

Wu, Y-Ta.; Nitschc, J.M. Chemical Engineering Science 
1995, 50, 1467-1487. 

Hanusse, P.; Blanche, A. J. Chem. Phys. 1981, 74, 6148. 



o 



o 



(a) 



o 



o 
o 



° o ° o 



o 



o 



o 



o 



o 



o o 




(b) 



FIG. 1. Various situations which axe simulated are shown. 
The three figures schematically depict various types of ini- 
tial conditions from which simulation is started, (a) The up- 
per most graph shows a situation where particles react when 
nearest-neighbors only. The reaction range is very short and 
particles come rarely in contact, (c) The lowest figure shows 
a situation of dense packing with a large reaction range. It 
corresponds to situation of high packing which occurs in a 
cell environment. It is unrealistic that particles can penetrate 
into each other but we consider this case nevertheless since it 
is simpler to model, (b) The middle graph is midway between 
two extremes. 
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FIG. 2. Result of Monte Carlo simulations in Id for type 

(i) initial condition. A very large system is simulated on 
a lattice with L = 10^ sites. Also, the initial number of 
particles Nq^a = Nq^b = 5000 is very large. Simulation 
starts from the largest possible density ntot(O) = 1 parti- 
cle/site. A and B particles have the same diffusion constant 
Da = Db = ls~^. Asymptotically, the number of particles 
approaches zero. There are three distinct regimes present; 
(a) of the mean field decay (— oo < logioit) < —2), (b) 
plateau where particle concentration does not change much 
(—2 < logio{t) < 2), (c) power law decay (2 < logioit) < 5), 
and (d) exponential decay at the end 5 < logioit) < oo. The 
indicated ranges are given roughly just to guide the eye. They 
also depend on which a is used in simulation. 
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FIG. 3. Study of the exponential regime where srnaU num- 
ber of particles is present on the lattice for type (i) initial con- 
dition with No, A = No,B = 5. The number of lattice sites is 
L — 1000. All other parameters are same as in figure 1. Each 
curve is obtained as average over 1000-3000 runs. Asymp- 
totically, the number of particles approaches zero. Panel (a) 
shows log-n versus log-t plot to trace down power law decay 
(should appear as a straight line). There is no power law de- 
cay. Also, small and largo t region are resolved better. Panel 
(b) shows log-n versus t plot to indicate exponential decay 
(corresponds to straight lines) . The particle density vanishes 
exponentially n ~ exp(— Kt) where k depends on a since the 
slopes for all curves are different. There is a value a = L/2 
when K becomes infinite (particles can not escape each other). 
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FIG. 4. Simulation for type (ii) initial condition. All pa- 
rameters as in the figure 2. The only difference from figure 2 
is in the initial condition. No. a and Nq,b vary randomly with 
constraint that A'o.a -|- No,b is fixed and equals 10. Asymp- 
totically, number of particles does not approach zero. 



-2 

-2.5 
-3 
-3.5 
o „4 

p 

-4.5 
-5 
-5.5 



(a) 



— Sim 

■ ■ ■ theory 
■ - theory+interp 




J i_ 



1 2 

iog,„(t) 




FIG. 5. Comparison of theory and experiment (Monte 
Carlo simulation) for type (i) initial condition. Simula- 
tion data are taken from figure 2. Panel (a): theory (dot- 
ted line, Eq. 15) predicts limt_>cx3 n{t) / while in reality 
n{t = oo) = 0. Reasons for this discrepancy are given in 
the text. Panel (b): by using interpolation formula (19) one 
obtains dashed curve. Agreement with simulation gets better. 
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FIG. 6. Comparison of theory and experiment (Monte 
Carlo simulation) for type (ii) initial condition. Both theory 

and simulation give n(oo) 7^ 0. Theory (dotted line, calcu- 
lated with Eq. 15) predicts faster annihilation of particles. 
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